%% 清空环境
clear; close all; clc;

%% ========== 参数与初值 (对应 Fig.8) ==========
a0 = -0.895;   % 耦合强度
a1 = 0.83;
a2 = 0.14;
b0 = 1.3;
b1 = 0.1;
b2 = 1;
d1 = 0;
d2 = 0;

x0 = 0.1;
y0 = 0.1;
z0 = 0.1;

N = 10336000;  % 迭代次数

%% ========== 预分配数组 ==========
x = zeros(1, N);
y = zeros(1, N);
z = zeros(1, N);

x(1) = x0;
y(1) = y0;
z(1) = z0;

%% ========== 迭代计算 ==========
for n = 1 : N-1
    x(n+1) = x(n) + y(n);
    y(n+1) = sin( y(n)*a0*sin(a1 + a2*x(n))*y(n) ) + b2*sin(z(n)) + d1;
    z(n+1) = b0*y(n) + b1*z(n) + d2;
end

%% ========== 绘制 Fig8 —— 时间序列 ==========
% 这里直接绘制 x 随迭代次数 t 的变化（红色点）
figure;
t = 1:N;
plot(t, x, 'r.', 'MarkerSize', 1);
xlabel('t');
ylabel('x_n');
title('Fig.8 时间序列图');
grid on;
hold on;

%% ========== 对 x 值进行聚类，确定 vortex 水平分界线 ==========
% 为加快计算，对 x 值进行降采样
decimation = 100;
x_sample = x(1:decimation:end);
% 设定聚类数 k = 6（根据论文描述有 6 种 vortex 状态）
k = 5;
% 使用 kmeans 聚类
[idx, centers] = kmeans(x_sample', k, 'Replicates', 10);
% 对聚类中心排序
centers = sort(centers);
% 计算相邻中心的中点，作为划分界线
boundaries = (centers(1:end-1) + centers(2:end)) / 2;

% 绘制水平分割线
yl = ylim;
for i = 1:length(boundaries)
    yline(boundaries(i), 'k--', 'LineWidth', 1);
end
hold off;

